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ABSTRACT 

Many countries are suffering from severe air pollution. Under- 
standing how different air pollutants accumulate and propagate is 
critical to making relevant public policies. In this paper, we use 
urban big data (air quality data and meteorological data) to iden- 
tify the spatiotemporal (ST) causal pathways for air pollutants. 
This problem is challenging because: (1) there are numerous noisy 
and low-pollution periods in the raw air quality data, which may 
lead to unreliable causality analysis; (2) for large-scale data in the 
ST space, the computational complexity of constructing a causal 
structure is very high; and (3) the ST causal pathways are com- 
plex due to the interactions of multiple pollutants and the influ- 
ence of environmental factors. Therefore, we present p-Causality , 
a novel pattern-aided causality analysis approach that combines the 
strengths of pattern mining and Bayesian learning to efficiently and 
faithfully identify the ST causal pathways. First, Pattern mining 
helps suppress the noise by capturing frequent evolving patterns 
(FEPs) of each monitoring sensor, and greatly reduce the complex- 
ity by selecting the pattern-matched sensors as “causers”. Then, 
Bayesian learning carefully encodes the local and ST causal re- 
lations with a Gaussian Bayesian network (GBN)-based graphical 
model, which also integrates environmental influences to minimize 
biases in the final results. We evaluate our approach with three real- 
world data sets containing 982 air quality sensors, in three regions 
of China from 01-Jun-2013 to 19-Dec-2015. Results show that 
our approach outperforms the traditional causal structure learning 
methods in time efficiency, inference accuracy and interpretability. 
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1. INTRODUCTION 

Recent years have witnessed the air pollution problem becoming 
a severe environmental and societal issue around the world. For 
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example, in 2015, the average concentration of PM2.5 in Beijing 
is greater than 150, classified as hazardous to human health by the 
World Health Organization, on more than 46 days. On Dec 7th 
2015, the Chinese government issues the first red alert because of 
the extremely heavy air pollution, leading to suspended schools, 
closed construction sites, and traffic restrictions. Though many 
ways have been deployed to reduce the air pollution, the severe 
air pollution in Beijing has not been significantly alleviated. 

Identifying the causalities has become an urgent problem for 
mitigating the air pollution and suggesting relevant public policy 
making. Previous research on the air pollution cause identifica- 
tion mostly relies on chemical receptor [1] or dispersion models 
[2]. However, these approaches often involve domain- specific data 
collection which is labor-intensive, or require theoretical assump- 
tions that real-world data may not guarantee. Recently, with the 
increasingly available air quality data collected by versatile sensors 
deployed in different regions, and pubic meteorological data, it is 
possible to analyze the causality of air pollution through a data- 
driven approach. 

The goal of our research is to learn the spatiotemporal (ST) causal 
pathways among different pollutants, by mining the dependencies 
among air pollutants under different environmental influences. Fig- 
ure 1 shows two example causal pathways for PM10 in Beijing. Let 
us first consider the pathway in Figure 1(a). When the wind speed is 
less than 5 m/s, the high concentration of PM10 in Beijing is mainly 
caused by SO 2 in Zhangjiakou and PM2.5 in Baoding. In contrast, 
as shown in Figure 1(b), when the wind speed is larger than 5m/s, 
PM10 in Beijing is mainly due to PM2.5 in Zhangjiakou and NO 2 
in Chengde. Based on this example, we can see the spatiotemporal 
(ST) causal pathways should reflect the following two aspects: 1) 
the structural dependency , which indicates the reactions and prop- 
agations of multiple pollutants in the ST space; and 2) the global 
confounder , which denotes how different environmental conditions 
could lead to different causal pathways. 




(a) Causal pathways (wind < 5m/s) (b) Causal pathways (wind > 5m/s) 

Figure 1: An illustration of identifying causal pathways. 

However, identifying the ST causal pathways from big air qual- 
ity and meteorological data is not trivial because of the following 
challenges. First , not all air pollution data are useful for causal- 


ity analysis. In the raw sensor-collected air quality data, there are 
numerous uninteresting fluctuations and noisy variations. Includ- 
ing such data into the causality analysis process is expected to lead 
to unreliable conclusions. Second , the sheer size of the air quality 
makes the causality analysis difficult. In most air quality moni- 
toring applications, thousands of sensors are deployed at different 
locations to record the air quality hourly for years. Discovering 
the ST causal relationships from such a large scale is challenging. 
Third , air pollution causal pathways are complex in nature. The 
air polluting process typically involves multiple types of pollutants 
that are mutually interacting, and is subject to local reactions, ST 
propagations and confounding factors, such as wind and humidity. 

Existing data mining techniques for learning the causal path- 
ways can be categorized into two classes: pattern-based [3] [4] and 
Bayesian-based [5] [6]. Pattern-based approaches aim to extract fre- 
quently occurring phenomena from historical data by applying pat- 
tern mining techniques; while Bayesian-based techniques use di- 
rected acyclic graphs (DAGs) to encode the causality and then learn 
the probabilistic dependencies from historical data. Though in- 
spiring results have been obtained by pattern-based and Bayesian- 
based techniques, both approaches have their merits and down- 
sides. Pattern-based approaches can fast extract a set of patterns 
(e.g., frequent patterns, contrast patterns) from historical air qual- 
ity data. Such patterns can capture the intrinsic regularity present 
in historical air quality data. However, they only provide shallow 
understanding of the air polluting process, and there are usually a 
huge number of frequent patterns, which largely limits the usability 
of the pattern set. On the other hand, Bayesian-based approaches 
depict the causal dependencies between multiple air pollutants in a 
principled way. However, the performance of Bayesian-based mod- 
els is highly dependent on the quality of the training data. When 
there exist massive noise and data sparsity, as the case of the air 
quality data, the performance of the Bayesian-based models is lim- 
ited. Besides, Bayesian-based approaches are limited by high com- 
putational cost [7] and the impact of confounding [8]. 

We propose p-Causality , which combines pattern mining with 
Bayesian learning to unleash the strengths of both. We claim p- 
Causality is essential for ST causal pathway identification, with 
the contributions listed as below: 

• First, a pattern mining module is designed to extract the fre- 
quent evolving patterns (FEPs) for each sensor, which helps sup- 
press the noise and reduce the computational complexity in the 
subsequent causal structure learning. 

• Second, a Bayesian learning module is proposed to systemati- 
cally encode and learn the causal pathways for air pollutants in the 
ST space, based on the high-quality data selected in the previous 
module. We also integrate the environmental factors to minimize 
the biases in the final results. 

• Third, we have carefully evaluated our proposed approach on 
three real data sets with 2.5 years’ air quality and meteorological 
data collected from hundreds of cities in China. Our results show 
that the proposed approach is significantly better than the existing 
baseline methods in time efficiency, inference accuracy and inter- 
pretability. 

2. FRAMEWORK 

In this section, we first describe the problem of identifying spatio- 
temporal causal pathways for air pollutants, and then introduce the 
framework of p-Causality. 

Let S = {s i , S 2 , • . • , sn, • • • } be the location set of the air qual- 
ity monitoring sensors deployed in a geographical region. Each 
sensor is deployed at a location s n £ S to periodically measure the 
target condition around it. All sensors have synchronized measure- 


ments over the time domain T = {1, 2, . . . , T}, where each t G T 
is a timestamp. We also consider a set C = {ci , C2 , . . . , cm } of pol- 
lutants. Given c m £ C, s n £ <S, and t £ T (1 < ra < M, 1 < n < 

N, 1 < t < T), we use P Crn s n t to denote the measurement of pol- 
lutant Cm at location s n and timestamp t. In addition, we also have 
the meteorological data at timestamp t for the entire geographical 
region, denoted as E t , as a vector of environmental factors. Using 
the air pollutant measurements and meteorological data, we aim to 
identify faithful causal relationships among different pollutants at 
different locations. 

Figure 2 shows the framework of our proposed approach p-Causality. 
It consists of two main modules: pattern mining and Bayesian Net- 
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Figure 2: The framework of our approach. 

Pattern Mining Module: This module first extracts th c frequent 
evolving patterns (FEPs) for each sensor. The FEPs essentially 
capture the air quality changing behaviors that frequently appear on 
the target sensor. By mining all FEPs from the historical air qual- 
ity data, this module efficiently captures the regularity in raw data 
and largely reduces the noise (Section 3.1 and 3.2). Afterwards, we 
examine the pattern-based similarities between locations to select 
candidate causers for each target sensor. By comparing the FEPs 
occurring on different sensors, we can obtain a shallow understand- 
ing of the causal relationships between different sensors, which can 
be further utilized to simplify learning the causal structures (Sec- 
tion 3.3). 

Bayesian Learning Module: By using the matched timestamps 
of the extracted FEPs at different sensors, together with the se- 
lected candidate sensors in the pattern mining module, this mod- 
ule further trains high-quality causal pathways from the large-scale 
air quality and context data in an effective and scalable way. We 
first generate the initial causal pathways from the selected candi- 
date causers, taking into account both the local interactions of mul- 
tiple air pollutants and the ST propagations (Section 4.1). Then to 
minimize the impact of confounding (Section 4.2), we integrate the 
confounders (e.g., wind, humidity) into the a Gaussian Bayesian 
Network (GBN)-based graphical model. Last, we refine the param- 
eters and structures of the Bayesian network to generate the final 
causal pathways (Section 4.3). 

3. THE PATTERN MINING MODULE 
3.1 Frequent Evolving Pattern 

To capture the frequent evolving behaviors of each sensor, we 
define frequent evolving pattern (FEP), an adaption of the clas- 
sic sequential pattern concept [9]. As the sequential patterns are 
defined on transactional sequences, we first discretize the raw air 
quality data. Given a pollutant c m at sensor s n , the measurements 
of Cm at s n over the time domain T form a time series. We dis- 
cretize the time series as follows: (1) partition it by day to obtain a 
collection of daily time series, denoted as P Crn s n ; and (2) for each 






daily time series {{pi, ti), (p 2 , t<i {pi, ti)), map every real- 
value measure pi (1 < i < l) to a discrete level pi using symbolic 
approximation aggregation [10]. After discretization, we obtain a 
database of symbolic sequences, as defined in Definition 1. 

Definition 1 (Symbolic Pollution Database). Forpol- 
lutant Cm and sensor s n , the symbolic pollution database P C ms n 
is a collection of daily sequences. Each sequence d G P Cm s n has 
the form {(pi, ti), (p 2 , £ 2 (pi, ti)) where an element (pi, ti) 
means the pollution level of c m at sensor s n and time ti is pi. 

Given the database P CmSn , our goal is to find frequent evolving 
behaviors of s n regarding c m . Below, we introduce the concepts of 
evolving sequence and occurrence. 

Definition 2 (Evolving Sequence). A length-k evolving 

sequence T has the form T = pi p 2 > • • • pk, where 

(l)\/i > l,pi-i pi and (2) At is the maximum transition time 
between consecutive records. 

Definitions (Occurrence). Given a daily sequence d = 
{(pi,ti), (P 2 , £ 2 ), • • • , (phti)) and an evolving sequence T = pi 
p 2 • • • pk (k < l), T occurs in d ( denoted as T □ d) if there 
exist integers 1 < ji < j ‘2 < • • • < jk < l such that: (1) VI < 
i < k, pj i = pi; and (2)\/l < i < k — 1, 0< tj i+1 — tj i < At. 

For clarity, we denote an evolving sequence pi -—$• p 2 • • • 

Pfe as pi -A p 2 • • • pk when the context is clear. Now, we 
proceed to define support wad frequent evolving pattern. 

Definition 4 (Support). Given P CrnSn and an evolving se- 
quence T, the support ofT is the number of days that T occurs, i.e., 
Sup(T) = |{o|o e Pc m s n ATE o}|. 

Definitions (Frequent Evolving Pattern). Given a sup- 
port threshold a, an evolving sequence T is a frequent evolving 
pattern in database Pcms n if Sup(T) > a. 

3.2 The FEP Mining Algorithm 

Now we proceed to discuss how to mine all FEPs in any symbolic 
pollution database. It is closely related to the classic sequential pat- 
tern mining problem. However, recall that there are two constraints 
in the definition of FEP: (1) the consecutive symbols must be differ- 
ent; and (2) the time gap between consecutive records should be no 
greater than the temporal constraint At. A sequential pattern min- 
ing algorithm needs to be tailored to ensure these two constraints 
are satisfied. 

We adapt PrefixSpan [9] as it has proved to be one of the most ef- 
ficient sequential pattern mining algorithms. The basic idea of Pre- 
fixSpan is to use short patterns as the prefix to project the database 
and progressively grow the short patterns by searching for local fre- 
quent items. For a short pattern ft, the /^-projected database Vp in- 
cludes the postfix from the sequences that contain fi. Local frequent 
items in Vp are then identified and appended to (3 to form longer 
patterns. Such a process is repeated recursively until no more local 
frequent items exist. One can refer to [9] for more details. 

Given a sequence a and a frequent item p, when creating p- 
projected database, the standard PrefixSpan procedure generates 
one postfix based on the first occurrence of p in a. This strategy, 
unfortunately, can miss FEPs in our problem. 

Example 1. Let At = 60 and a = 3. In the database shown 
in Table 1, item pi is frequent. The pi -projected database gener- 


Table 1: An example symbolic pollution database. 


Day 

Daily sequence 

di 

((P2,0), (Pi, 10), (p2, 30), (p3, 40)) 

d 2 

< (pi , 0) , (p 2 , 30) , (pi , 360) , (p 2 , 400) , (p 3 , 420)) 

d 3 

( (p2 , 0) , (p3, 30)) 

d/i 

<(pi, 0), (pi, 120), (p 3 , 140), (pa, 150), (p 3 , 180)) 

d 5 

<(P2, 50), (p 2 , 80), (p 3 , 120), (pi, 210)) 


ated by PrefixSpan is: 

(1) d,/p % m((p 2 ,20),(p 3 ,30)) 

(2) d 2 /p x = ((p 2 ,30), (pi, 360), (p2, 400), (p 3 , 420)) 

(3) d 4 /p x = ((pi, 120), (pa, 140), (pa, 150), (ps, 180)) 

The elements satisfying t < 60 are (p 2 , 20), (p 3 , 30) and (p 2 , 30). 
No local item is frequent, hence pi cannot be grown any more. 

To overcome this, given a sequence a and a frequent item p, we 
generate a postfix for every occurrence of p. 

Example 2. Also for Example 1 , if we generate a postfix for 
every occurrence of pi, the pi-projected database is: 

(1) di/pi = ((p2, 20), (j5 3 , 30)) 

(2) d 2 /p x = <(p2,30),(pi,360),(p2,400),(p3,420)> 

(3) da/ pi = ((P2,40),(p 3 ,60)) 

(4) d 4 /pi = ((pi, 120), (ps, 140), (p 2 , 150), (ps, 180)) 

(5) d 4 j p\ = <(P3,20),(P2,30),(P3,60)> 

The items p 2 and p 3 are frequent and meanwhile satisfy the tempo- 
ral constraint, thus longer patterns pi p 2 and pi p 3 are 

found in the projected database. 

Using the above projection principle, the projected database in- 
cludes all postfixes to avoid missing patterns under the time con- 
straint. Algorithm 1 sketches our algorithm for mining FEPs. The 
procedure is similar to the standard PrefixSpan algorithm in [9], 
except that the aforementioned full projection principle is adopted, 
and the time constraint At is checked when searching for local fre- 
quent items. 


Algorithm 1: Mining frequent evolving patterns. 


Input: support threshold cr, temporal constraint At, symbolic 



pollution database P 

1 Procedure InitialProjection( P , cr, At) 

2 

L frequent items in V; 

3 

foreach item i in L do 

4 


S 

5 


foreach sequence o in P do 

6 



R postfixes for all occurrences of i in o; 

7 



S <- S U R; 

8 


PrefixSpan(z, i, 1 , S, At); 

9 Function PrefixSpan(a, i pr ev, l, At) 

10 

L frequent items in S\ a meeting time constraint At; 

11 

foreach item i in L do 

12 


if 

i f iprev then 

13 



a! append i to a; 

14 



Build S\ a f using full projection; 

15 



Output a!; 

16 

_ 

- 

PrefixSpania', i, l + 1, S \ a r , At); 


The output of Algorithm 1 is the set of all FEPs for the given 
database, along with the occurring timestamps for each FEP. As an 





Figure 3: An illustration of the pattern-matched timestamps. 
The blue dashed lines represents the PM2.5 time series in Bei- 
jing during a two-year period, and the red points denote the 
timestamps at which a certain FEP has occurred (5 = 0.1). 


example, Figure 3 shows the raw PM2.5 time series in Beijing dur- 
ing a two-year period. After mining FEPs on the symbolic pollution 
database, we mark the timestamps at which the FEPs occur. One 
can observe that, the FEPs can effectively capture the regularly ap- 
pearing evolvements of PM2.5 in Beijing. Because of the support 
threshold and the evolving constraint, infrequent sudden changes 
and uninteresting fluctuations are all suppressed. 

3.3 Finding Candidate Causers 

After discovering the FEPs, next step is leverage them to extract 
the candidate causers for each sensor. Consider two sensors s and 
s' , let us use TS(s) and TS(s') to denote the sets of pattern start- 
ing timestamps for s and s', respectively. Below, we introduce the 
pattern match relationship. 

Definition 6 (Pattern Match). Let t s > e TS(s') be a 
timestamp at which a pattern happens on s'. For a pattern starting 
timestamp t s E TS(s), we say t 8 r matches t s if 0 < t s — t s r < L, 
where L is a pre-specified time lag threshold. 

Informally, the pattern match relation states that when there is 
a pattern occurring on s', then within some time interval, there is 
another pattern happening on s. Naturally, if s' has a strong causal 
effect on s, then most timestamps in TS s > will be matched by TS s , 
and vice versa. Based on TS s and TS s /, we proceed to introduce 
match precision and match recall to quantify the correlation be- 
tween s and s'. 


Definition 7 (Match Precision). Given TS s and TS s ’, 
we define the matched timestamp set ofTS s / as M s > — {t s / \ t s > E 
TS s ' A3t s E TS s , match(t s , t s r) = True}. With M s / andTS s >, 
we define the precision of s' matching s as: 

P( s ,s') = \M s t\/\TS a f\ 

Definition 8 (Match Recall). Given TS s and TS s r, we 
define the matched timestamp set ofTS s as M s — {t s \t s E TS s A 
3t s r E TS s r , match(t s ,t 8 /) = True}. With M s andTS s , we 
define the recall of s' matching s as: 

R{s,s') = \M s \/\TS s \ 


Relying on the concepts of match precision and match recall, 
we compute the pattern-based correlation between s and s' as the 
geometric mean of precision and recall, namely 


Corr(s , s') 


2 x P(s, s') 

P(s , s') + R(s, s') ’ 


Now we are ready to describe the process of finding candidate 
causers for each sensor. Given the set of all sensors and their 
pattern- starting timestamps, our goal is to find the candidate causers 
for each sensor. Consider a target sensor s, we say another sensor 


s' is a candidate causer for s if s' satisfies two constraints: (1) the 
distance between s and s' is no larger than a distance threshold S 9 ; 
and (2) the pattern correlation between s and s' is no less than a 
correlation threshold S p . Given the pattern- starting timestamps that 
are ordered chronologically, the retrieval of the candidate causers 
can be easily done by sequentially scanning the two timestamp lists 
to find pattern-matched pairs. 

Figure 4 illustrates eight examples of selected candidate causers. 
For PM2.5 in Beijing, we reduce the number of candidate sensors 
to X = 4 ~ 7 from overall |«S| = 61 sensors in North China. 
Note that China is a country with monsoon climate, the candidate 
sensors show quite similar geo-locations in four seasons. We there- 
fore separate the training data into four groups based on seasons, to 
better diagnose causalities for the air pollutants in China. 



Figure 4: Candidate sensors for Beijing PM2.5 in four seasons. 
Star: PM2.5 in Beijing. Circles: pollutants at candidate sen- 
sors. 


4. THE BAYESIAN LEARNING MODULE 
4.1 Generating Initial Causal Pathways 

To diagnose the complex reactions and dispersions of different 
air pollutants, both the local and ST dependencies need to be mod- 
elled in a fair way. Based on the pattern-matched data and the can- 
didate causers extracted by the pattern mining module, this subsec- 
tion first introduces the representation of causal pathways in the ST 
space, and then elaborates how to generate initial causal pathways. 

Definition 9 (Gaussian Bayesian Network (GBN)). 
GBN is a special form of Bayesian network for probabilistic infer- 
ence with continuous Gaussian variables in a DAG, in which each 
variable is assumed as linear function of its parents [11]. 


Local 

Qsot = ( Pc„s,(l-l) } 

Q(sFs N )t = { Pc m s n (t-l) } 
m = 1,2,...M; 

l = 1,2,...L; 

n = 1,2, ...,N 

(a) Local and ST dependencies in a GBN (b) Notations 

Figure 5: GBN-based causal pathway representation and its 
notations. 

As shown in Figure 5, the ST causal relations of air pollutants 
are encoded in a GBN-based graphical model, to represent both 
local and ST dependencies. Here we choose GBN to model the 
causalities because: 1) GBN provides a simple way to represent 
the dependencies among multiple pollutants variables, both locally 
and in the ST space. 2) GBN models continuous variables, while 




general Bayesian network (BN) models discrete variables. Since 
quantifying the continuous urban data into discrete values will in- 
troduce errors and more uncertainties for causality analysis, we use 
GBN instead of BN in our model. 3) The characteristics of ur- 
ban data fit the GBN model well, since the distribution of 1-hour 
difference (current value minus the value 1 -hour ago) of air pollu- 
tants and meteorological data obey Gaussian distribution (verified 
by D' Ago stino — Pearson test [12] [13]). In the following sec- 
tions, normalized 1 -hour differences of time series data will be used 
as inputs for the model. 

Specifically, for the target pollutant c m at sensor so-th sensor 
and timestamp t, denoted as P Crn s 0 t, m E [ 1 , M], we capture the 
dependencies from both the local causal pollutants Q^°£ al and 
the ST causal pollutants Here refer to 

a 1 x MNL vector of M types of pollutants at N neighborhood 
sensors s\ ~ sn and previous L timestamps that most probably 
cause the target pollutant in the ST space, i.e. — 

{ p c rn s ri (t-i)},rn = 1, . . . , M; n = 1, ...,N;Z = 1,...,L. Sim- 
ilarly, Qs 0 ° t caZ is a 1 x ML vector of pollutants locally at so- For 
example, when we set L as 2, M =. 6 , Qsot al m ay take values of 
12 normalized 1 -hour difference time series data, i.e. Q^°f al — 
(2, -0.5, 0.8, 0.3, 1, -2, 2.2, 1, 1,0, -0.5, 0.2). 

The parents of P Crn s 0 t are denoted as PA(P CrnSo t) = Qf o ° t cai 0 
Q(si~s N )t’ where ® denotes the concatenation operator for two 
vectors. Based on the definition of GBN, the distribution of P Crn s 0 t 
conditioned on PA(P CrnS0 t) obeys Gaussian distribution: 

Pr(PcmSot — PcmSot\PA(P CrnSo t)) ~ A/’(/X Cm s 0 t + 

P‘n=0 Pi l = l Pi m=l^ , m{nh-\-l){Pc rn s ri {t — l) Mc m s n (£ — l) ) 5 ^( £ c m sot)) 

Pcmsot is the marginal mean for P CrnSot . E denotes the Co- 
variance operator. A = {a m (nL + l)}, (m = 1, . . . , M; n — 

0, 1, . . . , N; / = 1, . . . , L) is the coefficient for the linear regres- 
sion in GBN [11]: 

To minimize the uncertainty of P Crn s 0 t given its parents, we need 
to find N sensors si ~ sn from the ST space and the parameters 
A that minimizes the error: 

S(e CmSot ) = £(P CmSot ) - AE(PA(P CmSot ))- 1 A T (2) 

Generating the initial causal pathways requires locating N most 
influential sensors from overall \S\ sensors requires up to (^) tri- 
als. Yet given the candidate sensors selected by Section 3.3, we 
manage to search from a subset ( X < \S\) sensors with time ef- 
ficiency and scalability. We further propose a a Granger causality 
score GC score to generate initial causal pathways, which is defined 
as: 

GC SCO re(rn, so, s n ) — maXm' g [ i,m] { \match(t ( Cm jS0 ) , t (c m / ,s n ) ) I 

I^( £ c m so0l | ~ |^( £ e m s 0 t)2| 'i 
P( £ c mSo t)2|x^(0.05) 1 

( 3 ) 

where GC score is a x 2 -test score [14] for the predictive causal- 
ity, with higher score indicating more probable “Granger” causes 
from M pollutants at sensor s n to the target pollutant c m at sensor 
so [15] ( GCscore < 1 means none causality). For variables obey- 
ing Gaussian distribution, Granger causality is in the same form 
as conditional mutual information [16], which has been used suc- 
cessfully for constructing structures for Bayesian networks. Here 
|raatc/i(£( CmjS0 ), t( c r Sn )) \ is the number of matched timestamps 
of FEPs between two time series (pollutant c m / at sensor s n and 
pollutant Cm at sensor so, see Section 3.3). And E(e Crn s 0 t)i and 
S(e Cm s 0 t )2 correspond to the variances of the target pollutant PcmSQt 


conditioned on Qi°£ al and Q^° f cal ® Q?J t . 

4.2 Integrating Confounders 

Recall the example in Figure 1. A target pollutant is likely to 
have several different causal pathways under different environmen- 
tal conditions, which indicates the environmental factors (humidity, 
wind, etc.) can be treated as confounders to the causalities among 
pollutants. This subsection is designed for integrating the environ- 
mental factors into the GBN-based graphical model, to minimize 
the biases in causality analysis and guarantee the causal pathways 
are faithful for the government’s decision making. We first intro- 
duce the definition of confounder and then elaborate the integration. 

Definition 10 (Confounder). A confounder is defined as 
a third variable that simultaneously correlates with the cause and 
effect, e.g. gender K may affect the effect of recovery P given a 
medicine Q, as shown in Figure 6(a). Ignoring the confounders will 
lead to biased causality analysis. To guarantee an unbiased causal 
inference, the cause-and-effect is usually adjusted by averaging 
all the sub- classification cases of K [5], i.e. Pr(P\do(Q)) m 
xZ =1 P r (P\Q,k)Pr(k). 
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(a) Cause-and-effect with confounder 
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(b) An illustration of cause-and-effect with confounders (environmental factors) 
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Figure 6: The GBN-based graphical model, integrating con- 
founders to the causal pathway, and converting the model into 
a generative model 

For integrating environmental factors, denoted as E t = {E ^ , E [ 2 ^ , 
into the GBN-based causal pathways, one challenge is there can 
be too many sub-classifications of environmental statuses. For ex- 
ample, if there are 5 environmental factors and each factor has 4 
statuses, there will exist 4 5 = 1024 causal pathways for each sub- 
classification case. Thus directly integrating E t as confounders 
to the cause and effect will result in unreliable causality analysis 
due to very few sample data conditioned on each sub-classification 
case. Therefore, we introduce a discrete hidden confounding vari- 
able K , which determines the probabilities of different causal path- 
ways from Q t to Pt , as shown in Figure 6(b). The environmental 
factors E t are further integrated into K, where K = 1, 2, ...K. 

In this ways, the large number of sub-classification cases of con- 
founders will be greatly reduced to a small number K, indicating K 
clusters of the environmental factors E t . 

Based on Markov equivalence (DAGs which share the same joint 
probability distribution [17]), we can reverse the arrow E t — » K 
to K E t , as shown in the right part of Figure 6(b). K de- 
termines the distributions of P, Qt,E t , thus enabling us to learn 
the distribution of the graphical model from a generative process. 

To help us learn the hidden variable K , the generative process fur- 
ther introduces a hyper-parameter n (as shown in Figure 6(c)) that 


determines the distribution of K. Thus the graphical model can 
be understood as a mixture model under K clusters. We learn the 
parameters of the graphical model by maximizing the new log like- 
lihood: 


4.3.2 M-step 

The membership probability in .E-step can be used to calculate 
new parameter values tt™ 6 ™, A £ ew , B % ew . We first determine 
the most likely assignment tag of timestamp t to cluster k , i.e. 


LL gen = Y lt T% =1 ln(Pr(pt\qt:k)Pr(qt\k)Pr(et\k)Pr(k\ir)) 

(4) 

In determining the number of the hidden variable K, we do not 
consider too large a K since that will induce much complexity for 
causality analysis. Also too small a K may not characterize the 
information contained in the confounders (i.e. meteorology). We 
observe the 2-D PCA projections of meteorological data (as shown 
in Figure 7). In three regions, five clusters can characterize the data 
sufficiently well. Thus we choose K = 3 ~ 7 for learning. 



(a) North China (NC) (b) Yangtze River Delta (YRD) (c) Pearl River Delta (PRD) 


Figure 7: 2-D PCA projections of 5 clusters of meteorological 
data in NC, YRD and PRD. The original meteorological data 
contains five types, i.e. temperature (T), pressure (P), humid- 
ity (U), wind speed (WS), and wind direction (WD), with each 
region divided into 9 grids, thus 45-dimensional. 


Tag t = max ke [ ljK ]7v t k (7) 

By integrating the timestamps belonging to each cluster k, we 
can update A k ew by Equation 5. Then we update B k by: 

TLCW f v-\^ r T 1 viT 

P*B k — 

* ( 8 ) 

ssr = ^T=nt k (e t - MSrXe* - m£7) t 

In addition, we update 7r^ k w by: 

new Itk 

7V t k =^tm ljT (9) 

SR phase (line 19-24) utilizes the parameters provided by the EM 
learning phase, and re-select the top N neighborhood sensors based 
on the newly generated GC score for each cluster k. We present a 
training example (as shown in Figure 8 (a)) of learning the causal 
pathways for Beijing PM2.5 during Jan^Mar. After 13 training 
iterations of the EM learning phase and structure reconstruction, 
we finally obtain K = 5 causal structures under each cluster, with 
the log likelihood shown in Figure 8 (b). For the last iteration, we 
calculate the percentage of labeled timestamps belonging to each 
cluster k. 


4.3 Refining Causal Structures 

This subsection tries to refine the causal structures and obtain 
the final causal structures under K clusters. The refining process 
includes two phases in each iteration: 1) an EM learning (EML) 
phase to infer the parameters of the model, and 2 ) a structure recon- 
struction (SR) phase to re- select the top N neighborhood sensors 
based on the new learnt parameters and GC SC ore , as illustrated in 
Algorithm 2. 

EML (line 6-18) is an approximation method to learn the param- 
eters 7r, 7 , Afc, B k of the graphical model, by maximizing the log 
likelihood (Equation 4) of the observed data sets via an .E-step and 
a M-step. Here 7 r contains the hyper parameters which determine 
the distribution of K (T x K-dimensional). 7 are posterior proba- 
bilities for each monitoring record (T x K-dimensional). A*., Bk 
are parameters for measuring the dependencies among pollutants 
and meteorology (K-dimensional). Note that A*., Bk come in dif- 
ferent formats. Ak is the regression parameter for: 

PcmS 0 t = ^0 + (Qiot al © Q(si~a N )t)Ak + e CmS 0 t (5) 

and Bk = (gs fc , Ss fe ) includes the parameters for the multi- 
variate Gaussian distribution of environmental factors E t . In the 
.E-step, we calculate the expectation of log likelihood (Equation 6 ) 
with the current parameters, and the M-step re-computes the pa- 
rameters. 


4.3.1 E-step 

Given the parameters 7 r, K, N, A&, Bk , EM assumes the mem- 
bership probability of pt,qt , e t belonging to the k - th cluster as: 


7 tk = Pr(k\p t ,qt,e t ) = 


Pr(k)Pr(p t ,qt,e t \k) 


Pr(p t ,qt,e t ) 
ntkJ\f(pt\qt, A k )JV \e t \Bk) 
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Algorithm 2: Refining the causal structures for each target pol- 
lutant Cm at location so- 

Input: T, K, N, and raining data sets pt, qt, et,t E [1, T] 

Output: Refined causal structures for K clusters 

Initial neighborhood sensors s\ ~ sn based on top N GC score', 

repeat 

EML(Pt, Qt, Et, si ~ K) 

Log Jikelihood, n tk , 7 tk , A k , B k ; 

SR(A fc , si ~ Sn , K)— »• s[ ~ s' N , Q'; 
until Log_likeoihood converges'. 

Function EM_Learning(EML)(Pt, Qt, Et, si ~ sjst, K) 

repeat 

Initial Assign: K clusters via K-means(Et) 
foreach item t = 1 to T do 

foreach item k = 1 to K do 
|_ Update 7 T t k by Equation (9); 

foreach item k = 1 to K do 
|_ Update Ak , Bk by Equation (5), (8); 

foreach item t = 1 to T do 

foreach item k = 1 to K do 
|_ Update 7 tk by Equation (6); 

until Log likelihood converges’, 

return: Log -likelihood and i Ttk , 'Jtk » Ak , Bk ; 

Function Structure -Reconstruction SR)(Ak, s i ~ sjv> K) 
foreach item s n in All candidate sensors do 

Compute GC score (jn, so, s n ) for si ~ sn’, 

Rank GC score and re-select the top N neighborhood sensors 

s i ~ S N’ 

Update Q — >• Q' corresponding to s [ ~ s' N ; 
return: s [ ~ s' N ,Q'’, 


5. EXPERIMENTS 



Figure 8: An example of learning the causal pathway for Beijing PM2.5, Jan^Mar under K — 5 clusters 


We evaluate the empirical performance of our method in this sec- 
tion. All the experiments were conducted on a computer with In- 
tel Core i5 3.3Ghz CPU and 16GB memory. We use MATLAB 
for our Bayesian learning module, and the open-source MATLAB 
BNT toolbox [18] for baseline methods. 

5.1 Experimental Setup 

5.1.1 DataSets 

We use three data sets that contain the records of 6 air pollutants 
and 5 meteorological measurements: 

• North China (NC), with 61 cities, 544 air quality monitoring 
sensors and 404 meteorological sensors in North China. The lati- 
tude and longitude ranges are 34N-43N, 110E-123E. 

• Yangtze River Delta (YRD), with 49 cities, 330 air quality 
monitoring sensors and 48 meteorology sensors. The latitude and 
longitude ranges are 28N-35N, 115E-123E, respectively. 

• Pearl River Delta (PRD), with 18 cities, 124 air quality moni- 
toring sensors and 406 meteorology sensors. The latitude and lon- 
gitude ranges are 22N-25N, 110E-116E. 

The 6 air pollutants are PM2.5, PM10, NO2, CO, O3, SO2, and 
the 5 meteorological measurements are temperature (T), pressure 
(P), humidity (H), wind speed (WS), and wind direction (WD). 
The time span for all data sets is from 06/01/2013 to 19/12/2015. 

We separate each data set into four groups based on four seasons, 
and use the last 15 days in each season in year 2014 and 2015 
for testing, and the remaining data for model training. We divide 
each region into 3x3 grids and average the meteorology values 
within each grid to get the environmental factors E t for the cou- 
pled model. 

We conduct experiments at both city level (Section 5.2.2, 5.2.1, 
5.2.4) and sensor level (Section 5.2.3). The city-level experiments 
average value of the sensors in the city to form a pseudo sensor, and 
discover the pathways among all the cities in three data sets. The 
sensor-level experiments analyze the causal relationships among 
sensors in each data set. 

5.1.2 Baselines 

Since Bayesian-based methods have been well used to learn causal 
Bayesian structures [19], we choose the most commonly used BN 
structure learning approaches as baselines to compare with our method. 
To identify the dependencies among different pollutants, the base- 
lines are deployed to learn the causal structures for each target pol- 
lutant. 

1. MWST. Maximum Weighted Spanning Tree (MWST) generates 
an undirected tree structure based on the MWST algorithm [20]. 
Each time it connects one edge between two nodes with the max- 
imum mutual information. Furthermore, [21] proposed an inde- 
pendency test method to assign a direction to each edge in the tree 
structure. 


2. MCMC. Markov-chain Monte Carlo (MCMC) is a statistical 
method that also samples from the Directed Acyclic Graph (DAG) 
space [22]. The method maximizes the score from a set of simi- 
lar DAGs that add, delete, or reverse connections, and updates the 
structure in the next iteration. 

2. K2+PS. K2 is a widely used greedy method for Bayesian struc- 
ture learning, which selects at most N parents based on the K 2 
score [23] for each variable given the updating order of all the vari- 
ables. In our case, we use pattern search algorithm [24] to optimize 
the updating order, thus reducing the search space of casual path- 
ways of different pollutants. Note that the original K 2 score is 
defined for discrete variables. Here we use GC score instead for the 
continuous variables. 

5.1.3 Parameter Setting 

The parameters of p-Causality are set as follows. (1) the sup- 
port threshold cr; (2) the temporal constraint At; (3) the distance 
threshold S 9 for finding candidate causers; and (4) the correlation 
threshold S p for finding candidate causers; (5) the number of time 
lags L = 3; (6) and the number of pollutant categories M = 6. When 
finding causal pathways at city level, we set cr = 0.1, At = 1 hour, 
S 9 = 200 km, and S p — 0.5. At the station level, all the the param- 
eters are set the same except that 5 g — 15 km to impose a finger 
granularity for finding candidate causers. K and N are evaluated 
within the range K = 3 ~ 7, and N= 1 ~ 5. 

5.2 Experimental Results 

5.2.1 Inference Accuracy 

We first evaluate the effectiveness of our approach via the causal 
inference accuracy through the causal pathways at city level, which 
is a 1-hour prediction task based on our proposed GBN-based graph- 
ical model. Note this prediction task is not general for all the 
timestamps, it only predict the future 1-hour based on the extracted 
pattern-matched periods, indicating the causal inference for the fre- 
quent evolving behaviors. Specifically, we first infer the probabil- 
ity Pr(k) of the testing data belong to cluster k. Then, we use the 
structure and parameters from the trained causal pathways regard- 
ing to this cluster to estimate the future pollutant concentration by 
Eq. 10. 

0 = + PA(P CmS0 t)A k )Pr(k) (10) 

The accuracy is defined as E^f* ( Pff 0 t-P2 rn s u t) latest, where 
Pcms 0 t is the ground truth value and T tes t is the number of test 
cases. Table 2 shows the 1-hour prediction accuracy for our ap- 
proach p-Causality, p-Causality-n which is p-Causality without in- 
tegrating the environmental factors, and three baseline methods in 
three regions of China. The accuracy shown in Table 2 is the av- 
erage accuracy for four seasons of each data set. The p-Causality 
gets the highest accuracy. The highest inference accuracy for the 


three regions are marked with three different colors (orange-> NC, 
blue->PRD, and green -> YRD) given different parameters N = 
3, K = 5. We note N = 3, K = 5 provides the best performance 
for Region NC, while N = 2, K = 3 generates the best accu- 
racy for YRD and N = 2, K = 5 for PRD. The optimal number 
N = 2 or N = 3 for the neighbor sensors also suggest that the 
air pollution is mainly affected by the most influential sensors in 
the ST space. We also find the accuracy of Region PRD is 76.7%, 
lower than Regions NC and YRD (86% and 89.1%). This may be 
explained by the fact that PRD region in China has relatively low 
pollution, while our approach targets at capturing the changes of 
pollutants and identifying their causality, thus may not be quite ap- 
plicable for prediction tasks of low pollution. We also evaluate the 
1-hour prediction accuracy with three well-used time series model, 
i.e., auto-regression moving average (ARMA) model, linear regres- 
sion model (LR), and support vector machine for regression with a 
Gaussian radial basis function (rbf) kernel (represented as SYM-R). 
p- Causality demonstrate more than 10% higher inference accuracy 
compared with these time series models. 


Table 2: Accuracy of PM2.5 1-hour prediction vs. baselines. 
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Table 3: Computation time for training data sets at city-level. 


Time (s) 

m = 1 

m = 2 

m = 3 

m = 4 

m = 5 

m = 6 

Mine FEPs 

2.84 

2.70 

3.34 

3.14 

3.47 

3.37 

Select candidate causers 

44.60 

61.26 

77.42 

55.25 

63.60 

62.50 

Generate initial causal pathways 

12.48 

13.51 

14.29 

15.32 

12.21 

13.60 

Integrate confounders 







Refine causal structures 

3954.96 

4316.40 

3661.68 

3987.12 

4507.68 

4267.92 

p-Causality 

4014.88 

4393.87 

3756.73 

4060.83 

4586.96 

4347.39 

Only Bayesian learning module 

4533.47 

4908.83 

4286.98 

4590.59 

5127.16 

4874.74 

MWST 

1572.29 

1709.04 

1538.4 

1591.1 

1796.48 

1787.25 

MCMC 

43548.28 

47291.59 

42053.29 

48314.53 

50290.17 

52014.68 

K2 + PS 

27290.18 

29378.45 

26355.63 

28213.09 

32058.1 

33248.55 


5 . 2.2 Time efficiency 

We also compare the training time of p- Causality with base- 
line methods. Since our approach consists of both pattern min- 
ing and Bayesian learning modules, we present the averaged time 
consumption of training all the three data sets, for each step in the 
two modules. We also evaluate the overall time consumption of 
p-Causality and p-Causality without pattern mining module (with 
only Bayesian learning). Results show that our approach is very ef- 
ficient, with the second minimum computation time among all the 
methods. MWST consumes the minimal time, however, it does not 
generate satisfactory accuracy for prediction (as in Section 5.2.1). 


We thus consider our approach provides the best trade-off regarding 
to accuracy and time efficiency. 

5.2.3 Scalability 

Another superior characteristic of our approach is the scalability. 
We further identify the causal pathways for air pollutants at sen- 
sor level, which is more than ten times as much in the city-level 
analysis. Our approach provides linear scalability in time with 1.1 
hours training time at city level, and 10 hours at sensor level. Note 
that when the data becomes larger, the pattern mining module can 
greatly save the computation time for generating the causal path- 
ways. Our approach needs about 10 hours, which is half of the 
training time of the approach without pattern mining. Moreover, 
the computation time can be greatly reduced compared to baseline 
methods. MCMC even cannot compute such large data sets and 
K2+PS takes up to 14 days, as shown in 9(b). Meanwhile, the ac- 
curacy is guaranteed when extending city-level data to sensor-level 
data, as shown in Figure 9(a). 



(a) Accuracy (b) Time efficiency 


Figure 9: Accuracy and time efficiency at station level 

5.2.4 Case Study 

To analyze the causal pathways for air pollutants, we study two 
cases corresponding to PM2.5 in specific cities. First we analyze 
the causal pathways for PM2.5 in the spring of Beijing and in the 
winter of Shanghai, the period of which are considered as the most 
heavily polluted season. Then we analyze Beijing PM2.5 before 
and during the APEC period (1st ~ 14th, Nov, 2014) as a case 
study for human intervention in causal systems. 

1. Beijing and Shanghai. Figure 8 is a real example for the 
causal pathways for Beijing PM2.5 during Jan^Mar. We provide 
the probability for each causal pathway for each cluster, defined 
as the proportion of labeled timestamps that belong to each clus- 
ter. As shown in Figure 8(a), Cluster 5 takes a dominant proportion 
(43.5%) of time for Beijing PM2.5, indicating the causal pathway 
during Jan^Mar most probably come from southern sensors, i.e. 
Baoding, Cangzhou and Langfang. Actions can be taken to control 
these pollutants in these cities. We then present the causal path- 
ways for PM2.5 in Shanghai, during Oct^Dec, which statistically 
has the highest air pollution concentration. As shown in Figure 
10, for PM2.5 in Shanghai, the N = 2 neighborhood cities gener- 
ally come from the northwest and the southwest. Cluster 2 takes 
a relative higher proportion (38.4%) of time for Shanghai PM2.5, 
suggesting the pollutants may be dispersed from PM2.5 in Nantong 
and S02 in Zhoushan. 

2. Beijing during APEC period. Traditionally, causality is veri- 
fied via interventions in a causal system. For example, we can ver- 
ify the effect of a medicine by setting two groups of patients and 
only giving medicine to the treatment group. However, it is impos- 
sible to conduct intervention for air pollutants in the real environ- 
ment. APEC period is a good opportunity to verify the causality, 
since the Chinese government shuttered factories in NC, and im- 
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Figure 10: Visualization of final causal pathways for PM2.5 in 
Shanghai. 


plemented traffic bans in and around Beijing [25]. Therefore, we 
compare the causal pathways for PM2.5 in Beijing before and dur- 
ing the APEC period. To illustrate the propagation of pollutants 
along the causal pathway, we connect the one-hop pathway to 3- 
hop as shown in Figure ll(a)(b). The connection originates from 
the target pollutant, i.e., Beijing PM2.5, and connect its causal pol- 
lutants at neighbor cities. Then for each new connected pollutant, 
we repeat the same procedure for the next hop. The connection 
stops if in inference accuracy of one target pollutant based on its 
historical data is higher than based on the historical data of its ST 
“causers”, indicating the pollutant is more likely to be generated 
locally. Figure 11(a) shows Beijing’s PM2.5 is likely to be caused 
by N02 in Zhangjiakou, Chengde and Baoding, 2 weeks before 
APEC period. Further, for example, Baoding ’s N02 is mostly in- 
fluenced by S02 in Shijiazhuang, Xingtai, and Hengshui. Note that 
the causal pathways forms “cliques” in the southwestern and north- 
eastern cities to Beijing, which is identical to the locations of the 
major plants in NC shown in Figure 11(c). However, we notice 
that the causal pathway cannot be connected into 3 -hops during the 
APEC period, since each “causer” pollutant to Beijing PM2.5 (i.e. 
N02 in Chengde and Zhangjiakou, PM10 in Cangzhou) is more 
likely to be inferred by its own historical data over its ST “causers” 
in this period. This may suggest the PM2.5 in Beijing during the 
APEC period are mostly affected by pollutants locally and nearby. 
The 3 -hop causal pathways learnt by three baselines are quite sim- 
ilar, thus we only present the result learned by MWST in Figure 
11(d). Our approach has better interpretability. We summarize the 
discovery for Beijing’s air pollution as follows. 
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(a) Causal pathway for Beijing PM2.5, (b) Causal pathway for Beijing PM2.5, 
Our method, two weeks before APEC Our method, APEC period 



(c) Locations of major plants in the 
center of North China 


(d) Causal pathway for Beijing 
PM2.5, MWST, Jan~Mar 


Figure 11: The causal pathways for Beijing PM2.5 before (a) 
and during APEC period (b), compared with the locations of 
major plants in Hebei Province, China (c), and the causal path- 
ways learned by baseline method MWST (d). 


• Among all the cities within a region, a target pollutant can be 
mainly affected by only several cities in the ST space. The locations 
of most influential cities to a target pollutant demonstrate seasonal 
similarities. 

• The causal pathways for PM2.5 in Beijing may come in mul- 
tiple hops that form “cliques” in certain regions, such as southwest 
and northeast of Beijing, suggesting superposition or reaction of 
air pollutants in the corresponding area. While during the APEC 
period with low pollution level, we did not see multi-hop causal 
pathways, suggesting the PM2.5 are more likely to be generated 
locally or nearby within this period. 

6. RELATED WORK 

In recent years, air pollution analysis has drawn a lot of atten- 
tion from the data mining community [26]. [27] [28] propose data- 
driven approaches to infer and forecast fine-grained air quality us- 
ing heterogeneous urban data. [29] estimates the gas consumption 
and pollutants emission of vehicles, based on the vehicles’ GPS tra- 
jectories in the road network. Our paper differs from these works 
in that, we focus on understanding the underlying causal pathways 
of air pollution, rather than predicting air quality or estimating pol- 
lutant emission. 

Causal modelling has been an active research direction during 
the past few years. Existing works on modeling causality for time 
series data can be classified into two categories. The first category 
considers a pair of time series, and aims to quantify the strength of 
causal influence from one time series to another. Researchers have 
developed different measures for this purpose, such as transfer en- 
tropy [16], Granger’s causality [15][14]. The second category aim 
to extract graphical causal relations from multiple time series. [30] 
combine graphical techniques with the classic Granger causality, 
and propose a model to infer causality strengths for a large num- 
ber of time series variables. Pearl’s causality model [5] encodes 
the causal relationships in a DAG [19] for probabilistic inference, 
with various extensions that incorporating spatial correlation [4] 
and temporal dependency [31]. 

Our proposed approach p-Causality belongs to the category of 
using graphical model to detect causalities from multiple time se- 
ries. It differs from the above works in two aspects: (1) We com- 
bine pattern mining and Bayesian learning to make the causality 
analysis more efficient and robust to the noise present in the in- 
put data. (2) Besides the multi- variate time series data, we also 
consider the impact of confounding given different environmental 
factors for unbiased causality analysis. 

7. CONCLUSION 

In this paper, we identified the ST causal pathways for air pollu- 
tants using large-scale air quality data and meteorological data. We 
have proposed a novel causal pathway learning approach named p- 
Causality that tightly combines pattern mining and Bayesian learn- 
ing. Specifically, by extending existing sequential pattern mining 
techniques, p-Causality first extracts a set of FEPs for each sen- 
sor, which capture most regularities in the air polluting process, 
largely suppress data noise and reduce the complexity in the ST 
space. In the Bayesian learning module, p-Causality leverages the 
pattern-matched data to train a graphical structure, which carefully 
models multi-faceted causality and environmental factors. We per- 
formed extensive experiments on three real- word data sets. Exper- 
imental results demonstrate that the causal pathways detected by 
p-Causality are highly interpretable and meaningful. Moreover, it 
outperforms baseline methods in both efficiency and inference ac- 
curacy. For future work, we plan to apply this pattern-aided causal- 


ity analysis framework for other tasks in the ST space, such as traf- 
fic congestion analysis. 

8. REFERENCES 

[1] S. Lee, W. Liu, Y. Wang, A. G. Russell, and E. S. Edgerton, 
“Source apportionment of pm 2.5: Comparing pmf and cmb 
results for four ambient monitoring sites in the southeastern 
united states ''Atmospheric Environment , vol. 42, no. 18, pp. 
4126-4137, 2008. 

[2] A. Keats, E. Yee, and F.-S. Lien, “Bayesian inference for 
source determination with applications to a complex urban 
environment,” Atmospheric environment , vol. 41, no. 3, pp. 
465-479, 2007. 

[3] C. Zhang, Y. Zheng, X. Ma, and J. Han, “Assembler: 

Efficient discovery of spatial co-evolving patterns in massive 
geo-sensory data,” in KDD. ACM, 2015, pp. 1415-1424. 

[4] H. Nguyen, W. Liu, and F. Chen, “Discovering congestion 
propagation patterns in spatio-temporal traffic data,” in the 
4th International Workshop on Urban Computing 

( UrbComp ), 2015. 

[5] J. Pearl, “Causality: models, reasoning and inference,” 
Economet. Theor , vol. 19, pp. 675-685, 2003. 

[6] J. Y. Zhu, Y. Zheng, X. Yi, and V. O. Li, “A gaussian 
bayesian model to identify spatio-temporal causalities for air 
pollution based on urban big data,” in Computer 
Communications Workshops (INFOCOM WKSHPS), 2016 
IEEE Conference on. IEEE, 2016. 

[7] D. M. Chickering, “Learning bayesian networks is 
np-complete,” in Learning from data. Springer, 1996, pp. 
121-130. 

[8] W. Sun, P. Wang, D. Yin, J. Yang, and Y. Chang, “Causal 
inference via sparse additive models with application to 
online advertising,” in AAAI, 2015, pp. 297-303. 

[9] J. Pei, J. Han, B. Mortazavi-Asl, H. Pinto, Q. Chen, 

U. Dayal, and M. Hsu, “Prefixspan: Mining sequential 
patterns by prefix-projected growth,” in ICDE , 2001, pp. 
215-224. 

[10] J. Lin, E. Keogh, S. Lonardi, and B. Chiu, “A symbolic 
representation of time series, with implications for streaming 
algorithms,” in SIGMOD, 2003. 

[11] M. A. Gomez, P. M. Villegasa, H. Navarrob, and R. Susia, 
“Dealing with uncertainty in gaussian bayesian networks 
from a regression perspective,” on Probabilistic Graphical 
Models , p. 145, 2010. 

[12] R. D’Agostino and E. S. Pearson, “Tests for departure from 
normality, empirical results for the distributions of b2 and 
bl,” Biometrika , vol. 60, no. 3, pp. 613-622, 1973. 

[13] J. Y. Zhu, C. Sun, and V. O. Li, “Granger-causality-based air 
quality estimation with spatio-temporal (st) heterogeneous 
big data,” in Computer Communications Workshops 

( INFOCOM WKSHPS), 2015 IEEE Conference on. IEEE, 
2015, pp. 612-617. 

[14] K. Hlavackova-Schindler, M. Palus, M. Vejmelka, and 
J. Bhattacharya, “Causality detection based on 
information-theoretic approaches in time series analysis,” 
Physics Reports , vol. 441, no. 1, pp. 1-46, 2007. 

[15] C. W. Granger, “Investigating causal relations by 
econometric models and cross-spectral methods,” 
Econometrica: Journal of the Econometric Society , pp. 
424-438, 1969. 

[16] L. Barnett, A. B. Barrett, and A. K. Seth, “Granger causality 


and transfer entropy are equivalent for gaussian variables,” 
Physical Review Letters , 2009. 

[17] I. Flesch and P. J. Lucas, “Markov equivalence in bayesian 
networks,” in Advances in Probabilistic Graphical Models. 
Springer, 2007, pp. 3-38. 

[18] K. Murphy et al. , “The bayes net toolbox for matlab,” 
Computing science and statistics , vol. 33, no. 2, pp. 
1024-1034,2001. 

[19] D. Heckerman, D. Geiger, and D. M. Chickering, “Learning 
bayesian networks: The combination of knowledge and 
statistical data,” Machine Learning , vol. 20, no. 3, pp. 
197-243, 1995. 

[20] C. Chow and C. Liu, “Approximating discrete probability 
distributions with dependence trees,” Information Theory, 
IEEE Transactions on , 1968. 

[21] G. Rebane and J. Pearl, “The recovery of causal poly-trees 
from statistical data,” arXiv preprint arXiv: 1304.2736, 2013. 

[22] J. L. Beck and S.-K. Au, “Bayesian updating of structural 
models and reliability using markov chain monte carlo 
simulation,” Journal of Engineering Mechanics , vol. 128, 
no. 4, pp. 380-391, 2002. 

[23] G. F. Cooper and E. Herskovits, “A bayesian method for the 
induction of probabilistic networks from data,” Machine 
Learning , vol. 9, no. 4, pp. 309-347, 1992. 

[24] R. M. Lewis and V. Torczon, “A globally convergent 
augmented lagrangian pattern search algorithm for 
optimization with general constraints and simple bounds,” 
SIAM Journal on Optimization , vol. 12, no. 4, pp. 

1075-1089, 2002. 

[25] K. Huang, X. Zhang, and Y. Lin, “The “apec blue” 
phenomenon: Regional emission control effects observed 
from space,” Atmospheric Research , vol. 164, pp. 65-75, 
2015. 

[26] Y. Zheng, L. Capra, O. Wolf son, and H. Yang, “Urban 
computing: concepts, methodologies, and applications,” 
ACM Transactions on Intelligent Systems and Technology 
(TIST), vol. 5, no. 3, p. 38, 2014. 

[27] Y. Zheng, F. Liu, and H.-P. Hsieh, “U-air: When urban air 
quality inference meets big data,” in KDD. ACM, 2013, pp. 
1436-1444. 

[28] Y. Zheng, X. Yi, M. Li, R. Li, Z. Shan, E. Chang, and T. Li, 
“Forecasting fine-grained air quality based on big data,” in 
KDD, 2015. 

[29] J. Shang, Y. Zheng, W. Tong, E. Chang, and Y. Yu, “Inferring 
gas consumption and pollution emission of vehicles 
throughout a city,” in KDD, 2014. 

[30] A. Arnold, Y. Liu, and N. Abe, “Temporal causal modeling 
with graphical granger methods,” in KDD. ACM, 2007, pp. 
66-75. 

[31] K. P. Murphy, “Dynamic bayesian networks,” Probabilistic 
Graphical Models, M. Jordan, vol. 7, 2002. 



